###############################################################################
## Replication Data for                                               #########
## Hemesath, S. & Tepe, M. (2025). Motivated Reasoning about          #########
##   Artificial Intelligence in Public Policy. Comparative Evidence   #########
##   from Germany. Behavioural Public Policy.                         #########
###############################################################################

## Load Packages

library(tidyverse)
library(qualtRics)
library(powerjoin)
library(survey)
library(marginaleffects)
library(ggh4x)

rm(list=ls())

## Read Data

df <- readRDS("Hemesath_Tepe_2025_Motivated_Reasoning_BPP.RDS")

###########################
#### Empirical Analysis ###
###########################

#### Code Gaais

gaaispos <- c("gaais1","gaais2","gaais4","gaais5","gaais7","gaais11","gaais12","gaais13","gaais14","gaais16","gaais17","gaais18")
gaaisneg <- c("gaais3","gaais6","gaais8","gaais9","gaais10","gaais15","gaais19","gaais20")

## Test internal consistency for each sample and each separate scale

# Sample 1
psych::alpha(df[df$Sample==1,][c(gaaisneg)])
psych::alpha(df[df$Sample==1,][c(gaaispos)])

# Sample 2
psych::alpha(df[df$Sample==2,][c(gaaisneg)])
psych::alpha(df[df$Sample==2,][c(gaaispos)])

# All samples

psych::alpha(df[c(gaaisneg)])
psych::alpha(df[c(gaaispos)])

### Create scale from both dimensions:

# Transform to numeric
df[c(gaaispos,gaaisneg)] <- df[c(gaaispos,gaaisneg)] %>% sapply(.,as.numeric)

## Internal consistency
psych::alpha(df[c(gaaispos)])
psych::alpha(df[c(gaaisneg)])

## predict factor scores
psych::fa(df[c(gaaispos)],nfactors=1)
psych::fa(df[c(gaaisneg)],nfactors=1)

df$gaais.pos <- psych::fa(df[c(gaaispos)],nfactors=1)$scores
df$gaais.neg <- psych::fa(df[c(gaaisneg)],nfactors=1)$scores

pcor <- ggcorrplot::ggcorrplot(as.data.frame(round(cor(df[c("gaais.pos","gaais.neg","techgov_1")]),2)), hc.order = TRUE, type = "lower",
           lab = TRUE)+
  theme_classic()+
  scale_x_discrete(labels=c("GAAIS(Negative)","Technology Regulation"))+
  scale_y_discrete(labels=c("GAAIS(Positive)","GAAIS(Negative)"))+
  theme(axis.title = element_blank())
pcor
ggsave(pcor,filename="Correlation.tiff",dpi=500,width=5,height=5)


## Attitudes towards Regulation of AI

# Reverse Code Items that measure more liberal positions towards AI
df$techgov_2 <- 6-df$techgov_2
df$techgov_4 <- 6-df$techgov_4

## test for Correlation and factorability of Items

psych::corr.test(df[,c("techgov_1","techgov_2","techgov_3","techgov_4")])
psych::KMO(df[,c("techgov_1","techgov_2","techgov_3","techgov_4")])

## Retest by Sample

psych::corr.test(df[df$Sample==1,][,c("techgov_1","techgov_2","techgov_3","techgov_4")])
psych::KMO(df[df$Sample==1,][,c("techgov_1","techgov_2","techgov_3","techgov_4")])

psych::corr.test(df[df$Sample==2,][,c("techgov_1","techgov_2","techgov_3","techgov_4")])
psych::KMO(df[df$Sample==2,][,c("techgov_1","techgov_2","techgov_3","techgov_4")])

# Plot distribution for GAAIS Scale and techgov_1 by Sample

p1 <- cowplot::plot_grid(ggplot(df %>% group_by(Sample,techgov_1) %>%
                                  summarize(n=n()) %>%
                                  mutate(freq=n/sum(n)),aes(x=techgov_1))+
                           geom_col(aes(y = freq*100,fill=Sample,color=Sample),
                                    position=position_dodge2(width=0.5)) +
                           geom_text(aes(y=(freq*100)+1,label=paste(round(freq*100,2),"%")),
                                     position=position_dodge2(width=1),color="grey",
                                     size=3)+
                           facet_wrap(.~"Attitudes towards Regulation of AI",)+
                           labs(x="Artificial intelligence should be more strictly\nregulated by the government.",
                                y="Share of responses")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1 (M=",round(mean(df$techgov_1[df$Sample==1]),2),")"),
                                                       paste("Sample 2 (M=",round(mean(df$techgov_1[df$Sample==2]),2),")")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1 (M=",round(mean(df$techgov_1[df$Sample==1]),2),")"),
                                                      paste("Sample 2 (M=",round(mean(df$techgov_1[df$Sample==2]),2),")")))+
                           theme_classic()+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank()),
                         ggplot(df,aes(x=gaais.pos))+
                           geom_density(alpha=0.1,linewidth=1,aes(color=Sample,fill=Sample))+
                           theme_classic()+
                           facet_wrap(.~"GAAIS\n(Index, Positive dimension)",)+
                           labs(x="Factor Scores",
                                y="Density")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1 (M=",round(mean(df$gaais.pos[df$Sample==1]),2),")"),
                                                       paste("Sample 2 (M=",round(mean(df$gaais.pos[df$Sample==2]),2),")")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1 (M=",round(mean(df$gaais.pos[df$Sample==1]),2),")"),
                                                      paste("Sample 2 (M=",round(mean(df$gaais.pos[df$Sample==2]),2),")")))+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank()),
                         ggplot(df,aes(x=gaais.neg))+
                           geom_density(alpha=0.1,linewidth=1,aes(color=Sample,fill=Sample))+
                           theme_classic()+
                           facet_wrap(.~"GAAIS\n(Index, negative dimension)",)+
                           labs(x="Factor Scores",
                                y="Density")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1 (M=",round(mean(df$gaais.neg[df$Sample==1]),2),")"),
                                                       paste("Sample 2 (M=",round(mean(df$gaais.neg[df$Sample==2]),2),")")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1 (M=",round(mean(df$gaais.neg[df$Sample==1]),2),")"),
                                                      paste("Sample 2 (M=",round(mean(df$gaais.neg[df$Sample==2]),2),")")))+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank()),
                         ggplot()+
                           theme_minimal())


p1
ggsave(p1,filename="Figure4.tiff",dpi=400,width=10,height=10,bg="white")

#######################################
### Socio-Demographic Distributions ###
#######################################

df$education <- as.factor(df$education)
levels(df$education) <- c("(Fach-)Abitur\n(High School Degree)",
                          "(Fach-)Abitur\n(High School Degree)",
                          "Hauptschule\n(Lower Secondary Degree)",
                          "(Applied) University Degree",
                          "Student",
                          "Without Degree",
                          "Realschule\n(Intermediate Degree)",
                          "(Applied) University Degree")
df$education <- fct_relevel(df$education,"(Applied) University Degree",after=0)
df$education <- fct_relevel(df$education,"Realschule\n(Intermediate Degree)",after=2)
df$education <- fct_rev(df$education)


px <- cowplot::plot_grid(ggplot(df %>% group_by(Sample,gender) %>%
                                  summarize(n=n()) %>%
                                  mutate(freq=n/sum(n)),aes(x=gender))+
                           geom_col(aes(y = freq*100,fill=Sample,color=Sample),
                                    position=position_dodge2(width=0.5)) +
                           geom_text(aes(y=(freq*100)+1,label=paste(round(freq*100,2),"%")),
                                     position=position_dodge2(width=1),color="grey",
                                     size=3)+
                           facet_wrap(.~"Gender",)+
                           labs(x="",
                                y="Share of respondents")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1"),
                                                       paste("Sample 2")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1"),
                                                      paste("Sample 2")))+
                           scale_x_discrete(labels=c("Non-Binary","Female","Male"))+
                           theme_classic()+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank()),
                         ggplot(df,aes(x=age))+
                           geom_density(alpha=0.1,linewidth=1,aes(color=Sample,fill=Sample))+
                           theme_classic()+
                           facet_wrap(.~"Age",)+
                           labs(x="Age",
                                y="Density")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1 (M=",round(mean(df$age[df$Sample==1]),2),")"),
                                                       paste("Sample 2 (M=",round(mean(df$age[df$Sample==2]),2),")")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1 (M=",round(mean(df$age[df$Sample==1]),2),")"),
                                                      paste("Sample 2 (M=",round(mean(df$age[df$Sample==2]),2),")")))+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank()),
                         ggplot(df %>% group_by(Sample,education) %>%
                                  summarize(n=n()) %>%
                                  mutate(freq=n/sum(n)),aes(x=education))+
                           geom_col(aes(y = freq*100,fill=Sample,color=Sample),
                                    position=position_dodge2(width=0.5)) +
                           geom_text(aes(y=(freq*100)+5,label=paste(round(freq*100,2),"%")),
                                     position=position_dodge2(width=1),color="grey",
                                     size=3)+
                           facet_wrap(.~"Education",)+
                           labs(x="",
                                y="Share of respondents")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1"),
                                                       paste("Sample 2")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1"),
                                                      paste("Sample 2")))+
                           theme_classic()+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank())+
                           coord_flip()+
                           scale_x_discrete(limits=rev),
                         ggplot(df %>% group_by(Sample,state) %>%
                                  summarize(n=n()) %>%
                                  mutate(freq=n/sum(n)),aes(x=state))+
                           geom_col(aes(y = freq*100,fill=Sample,color=Sample),
                                    position=position_dodge2(width=0.5)) +
                           geom_text(aes(y=(freq*100)+3,label=paste(round(freq*100,2),"%")),
                                     position=position_dodge2(width=1),color="grey",
                                     size=3)+
                           facet_wrap(.~"State",)+
                           labs(x="",
                                y="Share of respondents")+
                           scale_color_manual(values=c("orange","navy"),
                                              labels=c(paste("Sample 1"),
                                                       paste("Sample 2")))+
                           scale_fill_manual(values=c("orange","navy"),
                                             labels=c(paste("Sample 1"),
                                                      paste("Sample 2")))+
                           theme_classic()+
                           theme(legend.position = "bottom",
                                 legend.title=element_blank(),
                                 legend.background = element_blank())+
                           coord_flip()+
                           scale_x_discrete(limits=rev)

                         )


px
ggsave(px,filename="Appendix_Figure1.tiff",dpi=400,width=10,height=10)


### Test for order-effect of responses (by Sample)
df$MRtask[df$Sample==1] <- 1
df$MRtask <- as.factor(df$MRtask)

m1 <- glm(MRcorrectchoice ~ MRchoiceorder*Sample,data=df,family=binomial())

ame1 <- marginaleffects::marginaleffects(m1,variables=c("MRchoiceorder"),
                                         by=c("Sample"))


p2 <- ggplot(ame1,aes(x=contrast,y=estimate))+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,color=Sample,
                      alpha=ifelse(conf.low<0&conf.high<0,"a","b")),
                  position=position_dodge2(width=0.5),
                  fatten=8)+
  geom_hline(yintercept=0)+
  theme_classic()+
  facet_wrap("Order Effect of Correct Response Item\n(by Sample)"~.,)+
  scale_color_manual(values=c("orange","navy"))+
  theme(panel.grid.major.y = element_line(),
        legend.position="bottom")+
  scale_alpha_manual(values=c(1,.25))+
  guides(alpha="none",shape=guide_legend(title="Task"))+
  labs(y="Average Marginal Effect",x="")

p2
ggsave(p2,filename="OrderEffect.tiff",dpi=400,width=7,height=5)


#####################################################################
#### Baseline Results: Treatment Effect by level of congeniality ####
#####################################################################

### Code congeniality GAAIS

df$congeniality <- as.numeric(ifelse(df$MRTable == "A" | df$MRTable == "C" ,
                                     scales::rescale(df$gaais.pos,c(0,1)),
                                     scales::rescale(df$gaais.neg,c(0,1))))

### Code congeniality Tech Regulation

df$congeniality2 <- as.numeric(ifelse(
  df$MRTable == "A" | df$MRTable == "C",
  scales::rescale(df$techgov_1*-1,c(0,1)),
  scales::rescale(df$techgov_1,c(0,1))))


### Data Preparation

df$MRTable <- as.factor(as.character(df$MRTable))
df$MRContext <- as.factor(as.character(df$MRContext))
df$education_num <- as.numeric(df$education)


df$MRtreatment <- as.factor(ifelse(df$MRTable == "C" | df$MRTable == "D","Treatment","Placebo"))
df$MRContext <- fct_relevel(df$MRContext,"PublicAdministration",after=0)

m2 <- svyglm(MRcorrectchoice ~ congeniality*MRtreatment*MRContext+gender+education_num+age,
             family=binomial(link="logit"),
             design=svydesign(data=df,id=~ResponseId))
m3 <- svyglm(MRcorrectchoice ~ congeniality2*MRtreatment*MRContext+gender+education_num+age,
             family=binomial(link="logit"),
             design=svydesign(data=df,id=~ResponseId))


### Average Treatment Effect

## Estimate ATE

ame2 <- marginaleffects::marginaleffects(m2,
                                         variables=c("MRtreatment"),
                                         by=c("congeniality","MRContext"),
                                         newdata=datagrid(congeniality=seq(0,1,by=0.1),
                                                          MRContext = levels(df$MRContext)),
                                         vcov = "HC3"
                                         )

ame3 <- marginaleffects::marginaleffects(m3,
                                         variables=c("MRtreatment"),
                                         by=c("congeniality2","MRContext"),
                                         newdata=datagrid(congeniality2=seq(0,1,by=0.1),
                                                          MRContext=levels(df$MRContext)),
                                         vcov="HC3"
                                         )




# Plot Results by Context

# 1) Allocating Parking Permits

p3 <- cowplot::plot_grid(
  ggplot(ame3[ame3$MRContext=="PublicAdministration",],aes(x=congeniality2,y=estimate))+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom",
          legend.title=element_blank(),
          axis.title = element_text(face="italic"))+
    geom_smooth(aes(color=contrast))+
    facet_wrap("Allocating parking permits\n(Sample 1)"~.,)+
    geom_ribbon(aes(fill=contrast,ymin=conf.low,ymax=conf.high),
                alpha=0.4)+
    scale_color_manual(values=c("navy"))+
    scale_fill_manual(values=c("navy"))+
    labs(x="Congeniality\n(Technology Regulation)",
         y= "\u03B4 Predicted Correct Answer")+
    scale_y_continuous(limits=c(-.4,.41),breaks=seq(-1,1,by=.1))+
    geom_hline(yintercept=0),
  ggplot(ame2[ame2$MRContext=="PublicAdministration",],aes(x=congeniality,y=estimate))+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom",
          legend.title=element_blank(),
          axis.title = element_text(face="italic"))+
    geom_smooth(aes(color=contrast))+
    facet_wrap("Allocating parking permits\n(Sample 1)"~.,)+
    geom_ribbon(aes(fill=contrast,ymin=conf.low,ymax=conf.high),
                alpha=0.4)+
    scale_color_manual(values=c("navy"))+
    scale_fill_manual(values=c("navy"))+
    scale_y_continuous(limits=c(-.4,.41),breaks=seq(-1,1,by=.1))+
    labs(x="Congeniality\n(GAAIS)",
         y= "\u03B4 Predicted Correct Answer")+
    geom_hline(yintercept=0))
p3

# 2) Safety of Self-Driving Cars

p4 <- cowplot::plot_grid(
  ggplot(ame3[ame3$MRContext=="cars",],aes(x=congeniality2,y=estimate))+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom",
          legend.title=element_blank(),
          axis.title = element_text(face="italic"))+
    geom_smooth(aes(color=contrast))+
    facet_wrap("Safety of Self-Driving Cars\n(Sample 2)"~.,)+
    geom_ribbon(aes(fill=contrast,ymin=conf.low,ymax=conf.high),
                alpha=0.4)+
    scale_color_manual(values=c("navy"))+
    scale_fill_manual(values=c("navy"))+
    labs(x="Congeniality\n(Technology Regulation)",
         y= "\u03B4 Predicted Correct Answer")+
    scale_y_continuous(limits=c(-.3,.25),breaks=seq(-1,1,by=.1))+
    geom_hline(yintercept=0),
  ggplot(ame2[ame2$MRContext=="cars",],aes(x=congeniality,y=estimate))+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom",
          legend.title=element_blank(),
          axis.title = element_text(face="italic"))+
    geom_smooth(aes(color=contrast))+
    facet_wrap("Safety of Self-Driving Cars\n(Sample 2)"~.,)+
    geom_ribbon(aes(fill=contrast,ymin=conf.low,ymax=conf.high),
                alpha=0.4)+
    scale_color_manual(values=c("navy"))+
    scale_fill_manual(values=c("navy"))+
    scale_y_continuous(limits=c(-.3,.25),breaks=seq(-1,1,by=.1))+
    labs(x="Congeniality\n(GAAIS)",
         y= "\u03B4 Predicted Correct Answer")+
    geom_hline(yintercept=0))
p4

# 3) Prediction of Recidivism

p5 <- cowplot::plot_grid(
  ggplot(ame3[ame3$MRContext=="law" ,],aes(x=congeniality2,y=estimate))+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom",
          legend.title=element_blank(),
          axis.title = element_text(face="italic"))+
    geom_smooth(aes(color=contrast))+
    facet_wrap("Prediction of Recidivism\n(Sample2)"~.,)+
    geom_ribbon(aes(fill=contrast,ymin=conf.low,ymax=conf.high),
                alpha=0.4)+
    scale_color_manual(values=c("navy"))+
    scale_fill_manual(values=c("navy"))+
    scale_y_continuous(limits=c(-.3,.25),breaks=seq(-1,1,by=.1))+
    labs(x="Congeniality\n(Technology Regulation)",
         y= "\u03B4 Predicted Correct Answer")+
    geom_hline(yintercept=0),
  ggplot(ame2[ame2$MRContext=="law",],aes(x=congeniality,y=estimate))+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom",
          legend.title=element_blank(),
          axis.title = element_text(face="italic"))+
    geom_smooth(aes(color=contrast))+
    facet_wrap("Prediction of Recidivism\n(Sample2)"~.,)+
    geom_ribbon(aes(fill=contrast,ymin=conf.low,ymax=conf.high),
                alpha=0.4)+
    scale_color_manual(values=c("navy"))+
    scale_fill_manual(values=c("navy"))+
    scale_y_continuous(limits=c(-.3,.25),breaks=seq(-1,1,by=.1))+
    labs(x="Congeniality\n(GAAIS)",
         y= "\u03B4 Predicted Correct Answer")+
    geom_hline(yintercept=0))
p5

### Predicted Effects
## Estimate Marginal Effects for Treatment and Placebo Groups (by Context)

me2 <- marginaleffects::predictions(m2,
                                    vcov="HC3",
                                    newdata = datagrid(congeniality=seq(0,1,by=0.1),
                                                       MRContext=levels(df$MRContext),
                                                       MRtreatment=levels(df$MRtreatment))
                                    )
levels(me2$MRContext) <- c("Public Administration",
                           "Safety of Self-Driving Cars",
                           "Prediction of Recidivism")

me3 <- marginaleffects::predictions(m3,
                                    vcov="HC3",
                                    newdata = datagrid(congeniality2=seq(0,1,by=0.1),
                                                       MRContext=levels(df$MRContext),
                                                       MRtreatment=levels(df$MRtreatment))
                                    )
levels(me3$MRContext) <- c("Public Administration",
                           "Safety of Self-Driving Cars",
                           "Prediction of Recidivism")


## Plot Results

p6 <- cowplot::plot_grid(
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality2, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.2),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality2,y=estimate),
                data=me3[me3$MRContext=="Public Administration",])+
    facet_wrap("Allocating parking permits\n(Sample 1)"~.,)+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality2),
                alpha=0.4,
                data=me3[me3$MRContext=="Public Administration",])+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(Technology Regulation)",
         y= "Predicted Correct Answer")+
    geom_hline(yintercept=0.5)+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1))+
    guides(fill=guide_legend(title="Treatment Group",ncol=2),
           color=guide_legend(title="Treatment Group",ncol=2),
           linetype=guide_legend(title="Treatment Group")),
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.1),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality,y=estimate),
              data=me2[me2$MRContext=="Public Administration",])+
    facet_wrap("Allocating parking permits\n(Sample 1)"~.,)+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality),
                alpha=0.4,
                data=me2[me2$MRContext=="Public Administration",])+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(GAAIS)",
         y= "Predicted Correct Answer")+
    guides(fill=guide_legend(title="Treatment Group"),
           color=guide_legend(title="Treatment Group"),
           linetype=guide_legend(title="Treatment Group"))+
    geom_hline(yintercept=0.5)
)
p6

p7 <- cowplot::plot_grid(
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality2, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.2),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality2,y=estimate),
                data=me3[me3$MRContext=="Safety of Self-Driving Cars",])+
    facet_wrap("Safety of Self-Driving Cars\n(Sample2)"~.,)+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality2),
                alpha=0.4,
                data=me3[me3$MRContext=="Safety of Self-Driving Cars",])+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(Technology Regulation)",
         y= "Predicted Correct Answer")+
    geom_hline(yintercept=0.5)+
    scale_y_continuous(limits=c(0,.85),breaks=seq(0,1,by=.1))+
    guides(fill=guide_legend(title="Treatment Group",ncol=2),
           color=guide_legend(title="Treatment Group",ncol=2),
           linetype=guide_legend(title="Treatment Group")),
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.1),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality,y=estimate),
                data=me2[me2$MRContext=="Safety of Self-Driving Cars",])+
    facet_wrap("Safety of Self-Driving Cars\n(Sample2)"~.,)+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality),
                alpha=0.4,
                data=me2[me2$MRContext=="Safety of Self-Driving Cars",])+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_y_continuous(limits=c(0,.85),breaks=seq(0,1,by=.1))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(GAAIS)",
         y= "Predicted Correct Answer")+
    guides(fill=guide_legend(title="Treatment Group"),
           color=guide_legend(title="Treatment Group"),
           linetype=guide_legend(title="Treatment Group"))+
    geom_hline(yintercept=0.5)
)
p7

"Prediction of Recidivism"
"Prediction of Recidivism\n(Sample2)"

p8 <- cowplot::plot_grid(
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality2, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.2),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality2,y=estimate),
                data=me3[me3$MRContext=="Prediction of Recidivism",])+
    facet_wrap("Prediction of Recidivism\n(Sample2)"~.,)+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality2),
                alpha=0.4,
                data=me3[me3$MRContext=="Prediction of Recidivism",])+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(Technology Regulation)",
         y= "Predicted Correct Answer")+
    geom_hline(yintercept=0.5)+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1))+
    guides(fill=guide_legend(title="Treatment Group",ncol=2),
           color=guide_legend(title="Treatment Group",ncol=2),
           linetype=guide_legend(title="Treatment Group")),
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.1),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality,y=estimate),
                data=me2[me2$MRContext=="Prediction of Recidivism",])+
    facet_wrap("Prediction of Recidivism\n(Sample2)"~.,)+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality),
                alpha=0.4,
                data=me2[me2$MRContext=="Prediction of Recidivism",])+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(GAAIS)",
         y= "Predicted Correct Answer")+
    guides(fill=guide_legend(title="Treatment Group"),
           color=guide_legend(title="Treatment Group"),
           linetype=guide_legend(title="Treatment Group"))+
    geom_hline(yintercept=0.5)
)
p8


ggsave(p6,filename="Figure5.tiff",width=8,height=5)
ggsave(p7,filename="Figure6.tiff",width=8,height=5)
ggsave(p8,filename="Figure7.tiff",width=8,height=5)


############################################
### Robustness Tests (Subgroup Binning) ####
############################################

df$congeniality_bin <- factor(cut_number(df$congeniality,3))
df$congeniality2_bin <- as.factor(cut_number(df$congeniality2,3))

m2b <- svyglm(MRcorrectchoice ~ gender+education_num+age+congeniality_bin*MRtreatment*MRContext,
              family=binomial(link="logit"),
              design=svydesign(data=df,id=~ResponseId))
m3b <- svyglm(MRcorrectchoice ~ gender+education_num+age+congeniality2_bin*MRtreatment*MRContext,
              family=binomial(link="logit"),
              design=svydesign(data=df,id=~ResponseId))

### Average Treatment Effect

## Estimate ATE

ame2b <- marginaleffects::marginaleffects(m2b,
                                          variables=c("MRtreatment"),
                                          by=c("congeniality_bin","MRContext"),
                                          newdata=datagrid(congeniality_bin=levels(df$congeniality_bin),
                                                           MRContext = levels(df$MRContext)),
                                          vcov = "HC3"
)


ame3b <- marginaleffects::marginaleffects(m3b,
                                          variables=c("MRtreatment"),
                                          by=c("congeniality2_bin","MRContext"),
                                          newdata=datagrid(congeniality2_bin=levels(df$congeniality2_bin),
                                                           MRContext=levels(df$MRContext)),
                                          vcov="HC3"
)

ameb <- bind_rows(ame2b,ame3b%>%rename(congeniality_bin = congeniality2_bin),.id="ID")%>%
  mutate(.,CongenialityMeasure = factor(ID,labels=c("GAAIS",
                                                    "Technology Regulation")))
p9 <- ggplot(ameb %>%
         mutate(MRContext = factor(MRContext,
                                   labels=c("Allocation of Parking Permits\n(Sample 1)",
                                            "Safety of Self-Driving Cars\n(Sample 2)",
                                            "Prediction of Recidivism\n(Sample 2)")),
                CongenialityMeasure = factor(CongenialityMeasure,
                                             labels=c(
                                               "Congeniality based on\nAttitudes towards AI (GAAIS)",
                                               "Congeniality based on\nPreferences for Technology Regulation"
                                             ))),
       aes(x=congeniality_bin,y=estimate))+
  geom_hline(yintercept=0)+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,color = MRContext),
                  position=position_dodge2(width=.5))+
  facet_grid(.~fct_rev(CongenialityMeasure),scales="free_x")+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2))+
  scale_color_manual(values=c("navy","orange","purple"))+
  labs(x="Congeniality (Binned)",
       y="\u03B4 Predicted Correct Answer")+
  guides(color=guide_legend(title="Scenario"))
p9

ggsave(p9,filename="Appendix_Figure5.tiff",dpi=555,width=8,height=5)

###########################
## Robustness Moderators ##
###########################


df$gender <- as.factor(df$gender)
df$age_bin <- as.factor(cut_number(df$age,3))
df$LR_bin <- factor(cut_number(df$leftright_1,3),labels=c("Left 33%","Center 33%","Right 33%"))
df$education_bin <- factor(df$education,labels=c("Low","Low","Low","Medium","High","High"))


### GAAIS

### Gender

me2c <- marginaleffects::predictions(
    svyglm(MRcorrectchoice ~ education_num+age+congeniality_bin*MRtreatment*gender*MRContext,
           family=binomial(link="logit"),
           design=svydesign(data=df[df$gender != "Divers",],id=~ResponseId)),
                               vcov="HC3",
                               newdata = datagrid(congeniality_bin=levels(df$congeniality_bin),
                                                  MRtreatment=levels(df$MRtreatment),
                                                  gender=c("Frau","Mann"),
                                                  MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
           "Safety of Self-Driving Cars\n(Sample 2)",
           "Prediction of Recidivism\n(Sample 2)")))

p2c <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = gender,color=MRtreatment,
                      x=congeniality_bin,y=estimate),
                  data=me2c %>%
                    mutate(congeniality_bin=as.numeric(as.character(
                      factor(congeniality_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality,y=estimate),alpha=0.2,linewidth=0.1,
              data=me2%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality),
              alpha=0.1,data=me2%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
              )+
  facet_wrap(.~MRContext+MRtreatment,ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(17,16),labels=c("Female","Male"))+
  labs(x="Congeniality based on Attitudes towards AI (GAAIS)\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p2c

#Education

me2d <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+congeniality_bin*MRtreatment*education_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df[df$education_bin != "Divers",],id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality_bin=levels(df$congeniality_bin),
                     MRtreatment=levels(df$MRtreatment),
                     education_bin=levels(df$education_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p2d <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = education_bin,color=MRtreatment,
                      x=congeniality_bin,y=estimate),
                  data=me2d %>%
                    mutate(congeniality_bin=as.numeric(as.character(
                      factor(congeniality_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality,y=estimate),alpha=0.2,linewidth=0.1,
              data=me2%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality),
              alpha=0.1,data=me2%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards AI (GAAIS)\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p2d

# Age

me2e <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+congeniality_bin*MRtreatment*age_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df[df$age_bin != "Divers",],id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality_bin=levels(df$congeniality_bin),
                     MRtreatment=levels(df$MRtreatment),
                     age_bin=levels(df$age_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p2e <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = age_bin,color=MRtreatment,
                      x=congeniality_bin,y=estimate),
                  data=me2e %>%
                    mutate(congeniality_bin=as.numeric(as.character(
                      factor(congeniality_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality,y=estimate),alpha=0.2,linewidth=0.1,
              data=me2%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality),
              alpha=0.1,data=me2%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards AI (GAAIS)\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p2e

# LR

me2f <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+education_num+congeniality_bin*MRtreatment*LR_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df,id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality_bin=levels(df$congeniality_bin),
                     MRtreatment=levels(df$MRtreatment),
                     LR_bin=levels(df$LR_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p2f <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = LR_bin,color=MRtreatment,
                      x=congeniality_bin,y=estimate),
                  data=me2f %>%
                    mutate(congeniality_bin=as.numeric(as.character(
                      factor(congeniality_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality,y=estimate),alpha=0.2,linewidth=0.1,
              data=me2%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality),
              alpha=0.1,data=me2%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards AI (GAAIS)\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p2f


############################
### Technology Regulation ##
############################


### Gender

me3c <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ education_num+age+congeniality2_bin*MRtreatment*gender*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df[df$gender != "Divers",],id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality2_bin=levels(df$congeniality2_bin),
                     MRtreatment=levels(df$MRtreatment),
                     gender=c("Frau","Mann"),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p3c <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = gender,color=MRtreatment,
                      x=congeniality2_bin,y=estimate),
                  data=me3c %>%
                    mutate(congeniality2_bin=as.numeric(as.character(
                      factor(congeniality2_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality2,y=estimate),alpha=0.2,linewidth=0.1,
              data=me3%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality2),
              alpha=0.1,data=me3%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(17,16),labels=c("Female","Male"))+
  labs(x="Congeniality based on Attitudes towards the Regulation of AI\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p3c

#Education

me3d <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+congeniality2_bin*MRtreatment*education_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df[df$education_bin != "Divers",],id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality2_bin=levels(df$congeniality2_bin),
                     MRtreatment=levels(df$MRtreatment),
                     education_bin=levels(df$education_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p3d <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = education_bin,color=MRtreatment,
                      x=congeniality2_bin,y=estimate),
                  data=me3d %>%
                    mutate(congeniality2_bin=as.numeric(as.character(
                      factor(congeniality2_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality2,y=estimate),alpha=0.2,linewidth=0.1,
              data=me3%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality2),
              alpha=0.1,data=me3%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards the Regulation of AI\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p3d

# Age

me3e <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+congeniality2_bin*MRtreatment*age_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df[df$age_bin != "Divers",],id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality2_bin=levels(df$congeniality2_bin),
                     MRtreatment=levels(df$MRtreatment),
                     age_bin=levels(df$age_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p3e <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = age_bin,color=MRtreatment,
                      x=congeniality2_bin,y=estimate),
                  data=me3e %>%
                    mutate(congeniality2_bin=as.numeric(as.character(
                      factor(congeniality2_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality2,y=estimate),alpha=0.2,linewidth=0.1,
              data=me3%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality2),
              alpha=0.1,data=me3%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards the Regulation of AI\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p3e

# LR

me3f <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+education_num+congeniality2_bin*MRtreatment*LR_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df,id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality2_bin=levels(df$congeniality2_bin),
                     MRtreatment=levels(df$MRtreatment),
                     LR_bin=levels(df$LR_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p3f <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = LR_bin,color=MRtreatment,
                      x=congeniality2_bin,y=estimate),
                  data=me3f %>%
                    mutate(congeniality2_bin=as.numeric(as.character(
                      factor(congeniality2_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality2,y=estimate),alpha=0.2,linewidth=0.1,
              data=me3%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality2),
              alpha=0.1,data=me3%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards the Regulation of AI\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p3f


ggsave(p3e,file="Appendix_Figure6a.tiff",dpi=400,width=9,height=6)
ggsave(p3c,file="Appendix_Figure7a.tiff",dpi=400,width=9,height=6)
ggsave(p3d,file="Appendix_Figure8a.tiff",dpi=400,width=9,height=6)
ggsave(p3f,file="Appendix_Figure9a.tiff",dpi=400,width=9,height=6)

ggsave(p2e,file="Appendix_Figure6b.tiff",dpi=400,width=9,height=6)
ggsave(p2c,file="Appendix_Figure7b.tiff",dpi=400,width=9,height=6)
ggsave(p2d,file="Appendix_Figure8b.tiff",dpi=400,width=9,height=6)
ggsave(p2f,file="Appendix_Figure9b.tiff",dpi=400,width=9,height=6)

#################################
### Robustness: Response Time ###
#################################

df <- df %>% group_by(Sample)%>%
  mutate(duration_bin = factor(cut_number(`Duration (in seconds)`/median(`Duration (in seconds)`),3),
                               labels=c("Low","Medium","High")))

me2g <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+education_num+LR_bin+congeniality_bin*MRtreatment*duration_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df,id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality_bin=levels(df$congeniality_bin),
                     MRtreatment=levels(df$MRtreatment),
                     duration_bin=levels(df$duration_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p2g <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = duration_bin,color=MRtreatment,
                      x=congeniality_bin,y=estimate),
                  data=me2g %>%
                    mutate(congeniality_bin=as.numeric(as.character(
                      factor(congeniality_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality,y=estimate),alpha=0.2,linewidth=0.1,
              data=me2%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality),
              alpha=0.1,data=me2%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards AI (GAAIS)\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p2g


me3g <- marginaleffects::predictions(
  svyglm(MRcorrectchoice ~ gender+age+education_num+LR_bin+congeniality2_bin*MRtreatment*duration_bin*MRContext,
         family=binomial(link="logit"),
         design=svydesign(data=df,id=~ResponseId)),
  vcov="HC3",
  newdata = datagrid(congeniality2_bin=levels(df$congeniality2_bin),
                     MRtreatment=levels(df$MRtreatment),
                     duration_bin=levels(df$duration_bin),
                     MRContext=levels(df$MRContext)))%>%
  mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                               "Safety of Self-Driving Cars\n(Sample 2)",
                                               "Prediction of Recidivism\n(Sample 2)")))

p3g <- ggplot()+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high,shape = duration_bin,color=MRtreatment,
                      x=congeniality2_bin,y=estimate),
                  data=me3g %>%
                    mutate(congeniality2_bin=as.numeric(as.character(
                      factor(congeniality2_bin,labels=c("0.2","0.5","0.8"))))),
                  position=position_dodge2(width=0.25)) +
  geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,
                  x=congeniality2,y=estimate),alpha=0.2,linewidth=0.1,
              data=me3%>%
                mutate(MRContext = factor(MRContext,labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                             "Safety of Self-Driving Cars\n(Sample 2)",
                                                             "Prediction of Recidivism\n(Sample 2)"))))+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,y=estimate,x=congeniality2),
              alpha=0.1,data=me3%>%
                mutate(MRContext = factor(MRContext,
                                          labels=c("Allocation of Parking Permits\n(Sample 1)",
                                                   "Safety of Self-Driving Cars\n(Sample 2)",
                                                   "Prediction of Recidivism\n(Sample 2)")))
  )+
  facet_wrap(.~MRContext+MRtreatment,scales="free_x",ncol=4)+
  theme_classic()+
  theme(legend.position="bottom",
        panel.grid.major = element_line(),
        panel.grid.minor=element_line(linetype=2),
        legend.title=element_blank())+
  scale_color_manual(values=c("darkorange","navy"))+
  scale_fill_manual(values=c("darkorange","navy"))+
  scale_shape_manual(values=c(15,16,17))+
  labs(x="Congeniality based on Attitudes towards the Regulation of AI\n(Binned Terciles)",
       y="\u03B4 Predicted Correct Answer")+
  guides(fill="none",color="none",linetype="none")
p3g


ggsave(p2g,file="Appendix_Figure3.tiff",dpi=400,width=9,height=6)
ggsave(p3g,file="Appendix_Figure2.tiff",dpi=400,width=9,height=6)



########################################
### Secondary: AI Decides Themselves ###
########################################


### Predicted Effects
## Estimate Marginal Effects for Treatment and Placebo Groups (by Context)


df$techgov_4 <- 6-df$techgov_4

df$congeniality3 <- as.numeric(ifelse(
  df$MRTable == "A" | df$MRTable == "C",
  scales::rescale(df$techgov_4*-1,c(0,1)),
  scales::rescale(df$techgov_4,c(0,1))))

m7 <- svyglm(MRcorrectchoice ~ congeniality3*MRtreatment*MRContext+gender+education_num+age,
             family=binomial(link="logit"),
             design=svydesign(data=df,id=~ResponseId))

me7 <- marginaleffects::predictions(m7,
                                    vcov="HC3",
                                    newdata = datagrid(congeniality3=seq(0,1,by=0.1),
                                                       MRContext=levels(df$MRContext),
                                                       MRtreatment=levels(df$MRtreatment))
)

levels(me7$MRContext) <- c("Public Administration",
                           "Safety of Self-Driving Cars",
                           "Prediction of Recidivism")

## Plot Results

p7 <- cowplot::plot_grid(
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="none")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality3,y=estimate),
                data=me7[me7$MRContext=="Public Administration",])+
  geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality3),
              alpha=0.4,
              data=me7[me7$MRContext=="Public Administration",])+
  geom_col(data = df %>%
             mutate(bin = cut(congeniality3, breaks = 10)) %>%
             count(bin) %>%
             mutate(bin=seq(0.1,1,by=0.2),
                    frequency = n / sum(n)),
           aes(x = bin, y = frequency),
           alpha=0.25,
           fill="black")+
    facet_wrap("Allocating parking permits\n(Sample 1)"~.,)+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="",
         y= "Predicted Correct Answer")+
    geom_hline(yintercept=0.5)+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1)),
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="none")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality3,y=estimate),
                data=me7[me7$MRContext=="Safety of Self-Driving Cars",])+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality3),
                alpha=0.4,
                data=me7[me7$MRContext=="Safety of Self-Driving Cars",])+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality3, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.2),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
    facet_wrap("Safety of Self-Driving Cars\n(Sample 2)"~.,)+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="",
         y= "")+
    geom_hline(yintercept=0.5)+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1)),
  ggplot()+
    theme_classic()+
    theme(panel.grid.major = element_line(),
          legend.position="bottom")+
    geom_smooth(aes(color=MRtreatment,linetype=MRtreatment,x=congeniality3,y=estimate),
                data=me7[me7$MRContext=="Prediction of Recidivism",])+
    geom_ribbon(aes(fill=MRtreatment,ymin=conf.low,ymax=conf.high,x=congeniality3),
                alpha=0.4,
                data=me7[me7$MRContext=="Prediction of Recidivism",])+
    geom_col(data = df %>%
               mutate(bin = cut(congeniality3, breaks = 10)) %>%
               count(bin) %>%
               mutate(bin=seq(0.1,1,by=0.2),
                      frequency = n / sum(n)),
             aes(x = bin, y = frequency),
             alpha=0.25,
             fill="black")+
    facet_wrap("Prediction of Recidivism\n(Sample 2)"~.,)+
    scale_color_manual(values=c("darkorange","navy"))+
    scale_fill_manual(values=c("darkorange","navy"))+
    scale_linetype_manual(values=c(2,1))+
    labs(x="Congeniality\n(Each responsible themselves judging AI)",
         y= "Predicted Correct Answer")+
    geom_hline(yintercept=0.5)+
    scale_y_continuous(limits=c(0,0.8),breaks=seq(0,1,by=.1))+
    guides(fill=guide_legend(title="Treatment Group",ncol=2),
           color=guide_legend(title="Treatment Group",ncol=2),
           linetype=guide_legend(title="Treatment Group"))
)

levels(df$gender) <- c(2,1,3)
df$gender <- as.numeric(as.character(df$gender))

p7

ggsave(p7,file="Appendix_Figure4.tiff",dpi=400,width=9,height=6,bg="white")


### Export Regression Results in Table

## Rerun regression models with renamed variables for easier comparison




m2 <- svyglm(MRcorrectchoice ~ gender+education_num+age+congeniality*MRtreatment*MRContext,
             family=binomial(link="logit"),
             design=svydesign(data=df,id=~ResponseId))
m3 <- svyglm(MRcorrectchoice ~ gender+education_num+age+congeniality*MRtreatment*MRContext,
             family=binomial(link="logit"),
             design=svydesign(data=df %>% mutate(congeniality = congeniality2),id=~ResponseId))

m2b <- svyglm(MRcorrectchoice ~ gender+education_num+age+congeniality_bin*MRtreatment*MRContext,
              family=binomial(link="logit"),
              design=svydesign(data=df %>%
                                 mutate(congeniality_bin = factor(congeniality_bin,
                                                                  labels=c("Low","Medium","High"))),id=~ResponseId))
m3b <- svyglm(MRcorrectchoice ~ gender+education_num+age+congeniality_bin*MRtreatment*MRContext,
              family=binomial(link="logit"),
              design=svydesign(data=df %>%
                                 mutate(congeniality_bin = factor(congeniality2_bin,
                                                                  labels=c("Low","Medium","High"))),id=~ResponseId))


stargazer::stargazer(m2,
                     m3,
                     type="html",
                     column.labels = c("GAAIS",
                                       "Regulation"),
                     out="Baseline.html",
                     single.row = T)


stargazer::stargazer(m2b,
                     m3b,
                     type="html",
                     column.labels = c("GAAIS",
                                       "Regulation"),
                     out="Binned.html",
                     single.row = T)






